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Abstract 

Effective pair interactions with a soft-repulsive component are a well-known feature of polymer 
solutions and colloidal suspensions, but they also provide a key to interpret the high-pressure 
behaviour of simple elements. We have computed the zero-temperature phase diagram of four 
different model potentials with various degrees of core softness. Among the reviewed crystal struc- 
tures, there are also a number of non-Bravais lattices, chosen among those observed in real systems. 
Some of these crystals are indeed found to be stable for the selected potentials. We recognize an 
apparently universal trend for unbounded potentials, going from high- to low-coordinated crystal 
phases and back upon increasing the pressure. Conversely, a bounded repulsion may lead to in- 
termittent appearance of compact structures with compression and no eventual settling down in 
a specific phase. In both cases, the fluid phase repeatedly reenters at intermediate pressures, as 
suggested by a cell-theory treatment of the solids. These findings are of relevance for soft matter 
in general, but they also offer fresh insight into the mechanisms subtended to solid polymorphism 
in elemental substances. 

PACS numbers: 61.50.Ks, 61.66.Bi, 62.50.-p, 64.70.K- 
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I. INTRODUCTION 



Soft-matter systems, like solutions of polymer chains or dispersions of colloidal particles, 
have been the subject of increasing interest in the last few decades owing to the possibility 
they offer of exploring new kinds of equilibrium phase behaviour which radically depart from 
the simple-fluid paradigm as exemplified by rare gases . In the reduced-Hamiltonian 
approach, all of the colloid degrees of freedom are integrated out, except for the centre- 
of-mass coordinates which are taken to interact via an effective potential energy U. A 
further simplification occurs upon restricting the form of U to a sum of two-body spherically- 
symmetric contributions, U = J2i<j u ( r ij) > r ij being the distance between particles % and j. 
It turns out that the u(r) potential is typically softer (i.e., characterized by a less steep short- 
range repulsion) than, say, the Lennard- Jones potential, if not even being finite for vanishing 
interparticle distance as a result of full interpenetrability of the mesoscopic particles. 

A well-studied case is that of star polymers, i.e., colloidal particles with a more or less 
hard inner core surrounded by a soft penetrable corona made of grafted polymer chains. For 
those distances where coronas overlap, the effective repulsion between two stars grows only 
logarithmically with reducing distance, leading to stabilization of low-coordinated crystals 
in a range of densities . Another much-studied system is self-avoiding polymers, which 
are instead fully penetrable and roughly described by a simple Gaussian repulsion [5[. At low 
temperature, the Gaussian-core model (GCM) can exist in two distinct crystal phases of the 
face-centred cubic (fee) and body-centred cubic (bec) type 6, 7, Moreover, the fluid phase 
regains stability on increasing pressure at constant temperature (reentrant melting) due to 
an overwhelmingly larger number, at same energy, of spatially disordered configurations over 
crystalline ones. 

Simple elements under extreme thermodynamic conditions provide another instance of 
a soft short-range repulsion between the system constituents [9|. In this case, the softness 
of the effective repulsion is ultimately a reflection, on a coarse-grained level, of a pressure- 
dependent atomic radius, as determined in turn by pressure-induced charge transfer between 
atomic orbitals. These effects are well understood in the case of alkali metals, whose elec- 
tronic structure at high pressures departs radically from nearly free-electron behaviour and 
their common low-pressure symmetric structure (bec) becomes unstable to pairing of the 
ions. While it could generally be argued that the use of a classical interatomic potential is 
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not permitted for most elemental liquids and solids, owing to the fact that electronic effects 
are important and are even strongly enhanced by the pressure, a classical framework for 
the study of the phase behaviour of simple substances is still possible, in the same spirit 
of the Born-Oppenheimer separation between nuclear and electronic degrees of freedom in 
molecules. Non-adiabatic effects are truly negligible for insulators and semiconductors but 
seemingly not so for metals where the absence of a gap between occupied and unoccupied 
states would make the adiabatic approximation not particularly well justified. However, the 
fraction of electrons that can be scattered by phonons is rather small even at room temper- 
ature since electrons that lie sufficiently far from the Fermi surface remain frozen in their 
states by Pauli's principle. Hence, the notion of an adiabatic potential as well as the very 

also for normal metals. This conclusion can be 
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same concept of a crystal remain meaningfu 
made rigorous by a theorem due to Migdal 

In a recent publication 9J, we have sketched the overall phase diagram of a system of point 
particles interacting through the exp-6 (Buckingham JjJ) potential, which is being used since 
lon g ti me as an effective description of rare gases and metals under extreme conditions 12|, 
13l . LL4J . In spite of such popularity, important features of the exp-6 phase behaviour had 
previously passed unnoticed. In particular, the reentrant-melting behaviour, similar to the 
GCM model, and the rich solid polymorphism, both a generic trait of elemental substances 
at high pressures |15 |. In the exp-6 system, reentrant melting and solid polymorphism are 
both manifestations of the statistical competition, as a function of pressure, between two 
distinct scales of nearest-neighbour distance [9( . These two length scales result of 

the radial dependence of the exp-6 repulsion, which shows a range of distances where the 
force strength diminishes with decreasing interparticle separation (indeed, a rather strong 
form of core softness). Whence two preferred values of the mean neighbour distance and 
the ensuing frustration of crystalline packing of the standard (fee or bec) type in a certain 
pressure range. 

In this paper, we shall be concerned with the relation between softness of the interparticle 
repulsion and polymorphism of the solid phase, with an emphasis on the occurrence of 
thermodynamically stable non-Bravais crystals. We consider four distinct model potentials 
that have already appeared in the literature, so as to cover a wide range of possibilities. 
Aside from the degree of softness, the chosen potentials differ as for having or not a hard 
core as well as for diverging or not at zero separation. In all cases, we provide an accurate 
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analysis of the zero-temperature (T = 0) phase diagram by examining a large number of 
potentially relevant structures, much larger than ever considered before (at least for these 
potentials). In particular, we provide an update of the T = behaviour of the exp-6 
potential as being reported in Ref. We know of only a few studies of a similar kind 
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e structures was carried out more systematically (by e.g. a genetic 



201 ]. but neither regarded the potentials hereby analysed. 
The rest of the paper is organized as follows: In Section 2, we introduce the potentials 
for which the T = phases are computed. Then, in Section 3, we give a brief description 
of the crystal structures being analysed and briefly outline the numerical technique that we 
use to optimize a specific structure for a given pressure. Results for the different potentials 
are presented in Section 4, with an attempt to find a common thread among them. Some 
further remarks and conclusions are given in Section 5. 

II. SOFT-CORE POTENTIALS 

Model potentials that describe the effective pair interaction between particles of simple 
atomic fluids have a short-range repulsive component whose steepness (i.e., absolute slope) 
diverges at zero separation as an effect of the increasing hindrance, implied by the exclusion 
principle, to making particles more and more close to each other. In Lennard- Jones and 
inverse-power fluids, the repulsive force (= —u'(r)) steadily increases, with an ever increasing 
rate, as particles get closer. By contrast, in soft-core potentials the repulsive force, or at 
least the rate with which the force increases, is not a monotonous function of r. 

The first soft-core potential we want to consider is the exp-6 potential 

, . I +oo , r <a M 

u E {r) = < r «i (2.1) 

\^{6exp[a(l-r)]_ a (,) 6 } , r >a M 

where e is the depth of the attractive well, a is the position of the well minimum, a > 7 is a 
parameter governing the steepness of the short-range repulsion, and <jm is the point where 
the function in Eq. fl2.ll) attains its maximum. This interaction was already studied in some 



detail for the case a = 11 in [9], where we also provided a sketch of the phase diagram at 
T = 0. However, in that paper only a relatively small number of crystal structures were 
scrutinized. We here extend that study and find other crystals with a smaller chemical 
potential in a certain pressure range. 



The exp-6 potential as well as the Gaussian repulsion are instances of strongly-soft repul- 
sions in that the force slope is positive in a range of distances. However, we can also devise 
soft-core potentials whose force always increases upon reducing the interparticle distance 
while the force steepness does not monotonously increase as well. This is e.g. the case of 
the potential 
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where a = 2.1. This potential was first introduced in Refs. [21|, |22J, and will thus be referred 

to as the Yoshida-Kamakura potential. A limited study of the T = phase diagram of this 

model was already presented in Ref. 23[]. The force — u' Y K(r) always increases with reducing 

r but it does so at a somewhat smaller rate in a range of distances around a where u'y K (r) 

develops a local minimum (see Fig. 1). 

We consider two other model potentials that have recently been studied in the literature. 

These are the smoothed hard-core plus repulsive-step (Fomin) potential 24] . 
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with ro = 1.35 a, and the interaction potential of compressed elastic (Hertzian) spheres 25], 

K -> ' ~ (2.4) 
, r > a 

The main reason for considering these potentials is that an accurate analysis of the phase 
diagram was performed for both. Reentrant melting and solid polymorphism are observed in 
both cases but, as we shall see, not all stable low-temperature phases were actually identified. 
Moreover, while up{r) is a strongly-soft potential that, at variance with the exp-6 potential, 
lacks a strictly hard core, Mj^(r) is a bounded potential that, unlike the Gaussian repulsion, 
does not fit our definition of soft-core potential. This may appear strange since, for bounded 
repulsions, two particles can sit on top of each other - a fact which by itself is indication of 
the intrinsic "softness" of the interaction. 

To highlight the soft nature of the Hertz potential, we look at the r-dependence of what 
we call the effective inverse-power exponent (EIPE) of a potential. The EIPE of u(r) is 
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261 ] . By imposing the equality of the functions and their derivatives at a 
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In Fig. 2, we report n(r) for the potentials (I2.1l) - (l2.4p (for the sake of clarity, we have made 
the discontinuity of ue{t) at om milder by replacing the hard core with an extremely steep 
power-law barrier). We see indeed that, for all potentials, n(r) decreases significantly for 
decreasing r in a range of r values, thus signalling core softening and pointing at the same 
time to the existence of similarities among the four models as for the low- to moderate- 
pressure phase behaviour. 

We now prove, by adapting an argument originally due to Stillinger {(|, that a system 
of particles interacting through either uyk{ t ) or up(r) is found in a stable fee crystal for 
sufficiently low densities and temperatures, as indeed suggested by the divergence of the 
EIPE at infinity. Let u(r) = eip(r/a) be a purely repulsive and monotonously decreasing 
potential. The Boltmann's factor for a pair of particles is 



with R = ai/)~ 1 (kBT/e), ks being Boltzmann's constant (note that the distance R diverges 
as the temperature T goes to zero). We now calculate the T = limit of b(xR), 

fix) = lim b(xR) = exp {- lim . (2.7) 

For 

{exp{— 6(R/a) 2 \n(R/a)} (Yoshida — Kamakura potential) 
(2.8) 
exp{— 7_R/cr} (Fomin potential) , 

one invariably has 

0+ , < x < 1 
f(x) = { l/ e , x = 1 (2.9) 
1~ , x > 1 

In other words, at very low temperature u(r) reduces practically to the hard-sphere potential 
and exactly so in the limit T — > 0. Correspondingly, R plays the role of an effective hard-core 
diameter. Since the values of the hard-sphere packing fraction at freezing and melting are 
well known |27J, the lines of fluid-fee coexistence are implicitly given by the equations 

| (pi?(T) 3 ) / = 0.492 and | (pR(T) 3 ) m = 0.543 , (2.10) 

which become rigorously valid in the limit T — > 0. 
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III. CRYSTAL LATTICES AND THEIR SCRUTINY 



At T = and fixed pressure, a (crystal) phase is thermodynamically stable if its chemical 
potential is smaller than that of any other phase at equal pressure. We are then faced with 
the problem of minimizing the enthalpy per particle as a function of pressure among all 
possible crystals. This is a formidable problem, since the possibilities are virtually infinite. 
Hence, we restrict our search of stable structures to a limited albeit large number of candi- 
dates including, aside from Bravais crystals, also a number of Bravais lattices with a basis 
(i.e., non-Bravais lattices) that have been demonstrated to be relevant for some soft material 
or simple substance under high pressures. We do not anyhow consider the possibility for 
phases made up of clusters, columns, or lamellae, which can appear if the potential exhibits 
two competing length scales of considerable difference (see e.g. Ref. [3]). We are fully aware 
of the limitations of our approach which, besides assuming a preselected set of structures, 



also undergoes increasing difficulties when managing lattices with many (five or more^ 
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321 ] as a tool for an automated search of optimal solutions inside 



the structure space would ensure a higher rate of success. 

However, provided the number of crystal lattices being examined is sufficiently large, 
direct optimization of structures from a finite set offers the advantage of weighing up the 
relative stability between the optimum solution for a given pressure and the metastable 
crystals that crowd around it. Moreover, features like the trend of particle coordination 
as a function of pressure or the characteristic distances at which nearest and next-nearest 
particles are preferentially located around a central particle can anyway be grasped by this 
method. 

The crystal lattices that we analyse can be divided in groups of increasing optimization 
difficulty, according to the number (from zero to five) of independent parameters (i.e., axial 
ratios, angles, and/or atomic-site parameters), besides the number density p, that need to 
be adjusted in order to minimize the enthalpy at fixed pressure. For non-Bravais lattices, 
these parameters are listed in Table 1, along with the number of inequivalent sites (NIS) of 
the lattice, defined as the maximum number of sites whose environments look different as 
for the population and/or radius of at least one coordination shell. For a Bravais lattice one 
has NIS=1, but there also exist non-Bravais lattices with this property (e.g. diamond). 
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The simplest close-packed structures, fee and hexagonal close packed (hep), are well 
known. In the [111] direction, the stacking sequence of triangular layers is ABC for fee and 
AB for hep (the bee lattice can also be built this way, by suitably spacing the layers one 
relative to the other in the same ABC sequence as for fee). More generally, one can consider 
so-called Barlow packings (i.e., the stacking variants of fee and hep packings) where each 
additional layer involves a binary choice for how to place it relative to the previous layer. 
Among them, we consider the double hep (dhep) structure (stacking sequence: ABAC), 
which is the structure of aLa, the triple hep (thep) structure (stacking sequence: ABCACB), 



and the 9R structure of 5Sm and aLi 



36| (stacking sequence: ABCBCACAB). Another 



simple structure with no parameters is the diamond lattice (a fee lattice with a basis of two 
atoms), which provides the low-pressure ground-state configuration of C, Si, Ge, and Sn. 
Structures with zero parameters are also the simple-cubic (sc) lattice, the one-species analog 
of the fluorite (CaF 2 ) lattice, the A15 lattice (i.e., the structure of /3W, whose conventional 
unit cell contains eight atoms), and the bccl2 lattice (i.e., the structure of the metastable 
Ga-II phase). 

The group of lattices with one free parameter includes, besides a few Bravais lattices 
(sh, st, bet, and trig - i.e., the simple rhombohedral or hRl lattice), also the /3Sn lattice 
(occurring also for Si, Ge, Rb, and Cs), the graphite lattice, as well as two cubic lattices 
obtained by suitably distorting a bee supercell (cI16-Li [stJ and BC8, both with a conven- 
tional unit cell containing eight atoms). hRl is obtained by stretching the sc lattice along a 
body diagonal. The straining parameter h is usually defined [33} in such a way that h = 
gives the sc lattice, h = —1/6 the bee lattice, and h = 1/3 the fee lattice. Both Hg and 
Li display a hRl solid phase. The conventional unit cell of /3Sn is tetragonal with 4 atoms 



38. 



39| while BC8 provides a 



inside. A cI16 phase has recently been identified also for Na 
low-pressure metastable phase for Si. The coordinates of the atoms inside the conventional 
unit cells of cI16-Li and BC8 are written in terms of a parameter x ranging from (giving 
back the bee lattice in the cI16-Li case) to 1/8. 

Aside from the Bravais lattices so, sfco, bco, and fco, crystal structures with two ad- 
justable parameters are also the one-species analog of wurtzite (ZnS) and the so-called 
Imma phase of Si (Si-XI, stable between 13 and 16 GPa). The Imma phase is a distortion 
of /?Sn. Moreover, we consider: the /3Np structure, which is the one-species analog of PbO 
(Ad, tP4); the aAs structure (A7, hR2), also relevant for Sb and Bi, which can also be viewed 
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as an orthorhombic lattice with b/a = a/3 and a 12-atom basis; the 7Se structure (A8, hP3), 
also observed in Te; the /?Mn structure (A13, cP20); and, finally, the C6 structure of Ti, Zr, 
and Hf (lo phase, hP3). 

There are two further Bravais lattices with three parameters (sm, sfcm) and one with five 
parameters (trie). We also consider another three-parameter lattice (A20, oC4, Cmcm space 
group), which provides the structure of all, (3G&, and 7H, a four-parameter lattice (All, 
0C8, Cmca space group), and two more five-parameter lattices (oC16-Cs and ST12). The 
All lattice is the structure of aGa (black P and B above 90 GPa have a similar structure). 
A 0C8 phase was first predicted for compressed Li by Neaton and Ashcroft 40j and later 
confirmed by Hanfland et al. |37j. Theoretical calculations by Christensen and Novikov 
have predicted the existence of a stable 0C8 phase also for Na above ~ 250 GPa 41] . The 
orthorhombic 0CI6 lattice provides the structure of Cs-V [3_5j and Rb-VI but it is also 
observed in Si (between 38 and 45 GPa) and in Ge. Its Cmca structure contains two types 
of atom, say Ai and A2, with Ai in planar arrangements separating A2 double layers. Within 
the planes, atoms form a dense packing of dimers. The Cmca structure of 0C8 is similar to 
0CI6 but the double layers of A2 atoms are absent. Of Cmca symmetry is also the structure 
of Ca-V [42j, which has the highest superconducting transition temperature (25 K) among 
all the elements. Finally, ST12 denotes a simple-tetragonal metastable phase of Si with a 
12-atom basis. 

Let us now briefly outline the procedure we follow in order to optimize a given crystal 
structure for assigned potential. Loosely speaking, we make a grid of points in parameter 
space on which we search for the minimum enthalpy at given pressure. We consider lattices 
with a number of sites between 8000 and 10000, with periodic boundary conditions. Lattice 
sums are extended up to a cutoff distance equal to half the shortest box length for orthogonal 
cells (half the shortest distance of centre to boundary for non-orthogonal cells). Then, a 
long-range correction is added to the energy by assuming a radial distribution function 
(RDF) of 1 for the distant sites. We start by performing a series of harsh minimizations on 
rather coarse meshes but then we refine the calculation by progressively reducing the mesh 
size until the lattice providing the minimum chemical potential is extracted out of a number 
of deeper valleys in parameter space (steepest-descent methods would not necessarily be of 
help in this case since one could easily get stuck in local minima). We feel satisfied when 
we obtain the density with three exact decimal digits and the other parameters with two to 
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three decimals. Occasionally, we resolve near degeneracies by going to larger lattices and 
better parameter refining. When more than one crystal structure is found to give practically 
the same minimum chemical potential, we check the identity of the subtended lattices by 
looking at the discrete RDF of each. 



IV. RESULTS AND DISCUSSION 



The T = phase diagram of model systems interacting through the soft-core potentials 
f l2.ip - fl2.4|) is reported in Tables 2 to 5. Overall, we see the existence of a rich solid poly- 
morphism, with many exotic lattices providing stable phases at moderately high pressures. 
Obviously, we cannot exclude that there might exist other phases which overcome in stability 
some of those found (to be sure, one should perform a molecular-dynamics simulation of the 
system at fixed pressure, starting from the high-temperature fluid and cooling it down very 
slowly until the solid nucleates - but the conditions for observing ideal freezing are hardly 
achieved in a reasonable amount of computer time). In this respect, the case of Na is partic- 
ularly illuminating 39]: if the experiment does not lie, no theoretical study could ever have 
anticipated the bunch of crystal phases, some of which extremely complex, that Na exhibits 
for pressures above 110 GPa. Yet, we believe that our study gives some general teachings 
on the mechanisms underlying spatial ordering at high pressures which would remain valid 
even in case some of the phases that we find stable are actually metastable. 

As a preliminary test, we study the Lennard- Jones fluid at T = and confirm that it 
can only exist in two forms, either as a hep or as a fee crystal, the former being stable at 
lower pressures [43]. The transition between the two occurs approximately at Pa 3 /e = 800, 
though the chemical-potential difference is very small (< 0.001 e) for all pressures up to 
1000. Moreover, we find that apparently all Barlow packings are more stable than fee for 
Pa 3 /e < 800 while the metastable phase that is closest in stability to fee is cI16 (with 
x = 0.024 at P = 0, slightly diminishing on increasing pressure), as suggested by the 



simulation 



44|. We also verified that the only T = phases of the GCM are fee and bcc. 



Table 2 collects the T = phases of the exp-6 model with a — 11. In Fig. 3, the chemical 
potential fi of the same phases is plotted as a function of pressure P, chosen fee as reference 
(the units of length and energy are set to a and e, respectively). We do not show /x for all 
the scrutinized lattices since the difference in chemical potential between the stable and the 
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metastable phases is often very small on the scale of the picture, insomuch that it would 
have been difficult to keep track of all the curves. Examples of a strong competition between 
different phases are seen anyway in Fig. 3, near P = 15000 (hRl vs. sh) and near P = 35000 
(the challenge now being among oC8, hRl, sh, and wurtzite). 

It transpires from Table 2 that the coordination number z has a regular trend with 
pressure. Starting from 12 (fcc-I) at low pressure, z reduces upon compression down to a 
minimum of 2 (hRl). Then, it increases progressively with pressure until becoming, even- 
tually, 12 again (fcc-II). This behaviour can be rationalized as follows. We see from Table 
2 that the nearest-neighbour distance is om for all stable phases except for fcc-I and bcc-I. 
Whence the convenience, in order to minimize the energy at not too low pressure, that the 
number of first neighbours be as small as possible. Eventually, however, the PV term in the 
enthalpy takes over the energy and there will again be room for close-packed lattices. The 
reason why low-coordinated phases do not show up for systems with a Lennard- Jones-like 
potential is that, in these systems, the fee and bec crystals anyway manage to accomodate 
second and third neighbours at convenient distances. On the contrary, in systems interact- 
ing through soft-core potentials, fee and bec local orders are destabilized by the peculiar 
dependence of the interatomic force with distance, which, for high softness degrees, leads to 
the existence of two incommensurate length scales that heavily frustrate the too compact 
arrangements. 

Looking at Fig. 3, a few other comments are in order: 1) the BC8 phase, which was 
previously |9|] found to be stable for pressures between roughly 20000 and 30000, is actually 
less stable than /3Sn. 2) Interestingly, there are two distinct ranges of pressure where the 
exp-6 system exhibits oC8 order at T = 0. In fact, within one of these intervals there are 
phase transitions between different oC8 phases, as signalled by the jumps of parameters 
from one valley to another in parameter space. We tabulated up to four distinct oC8 phases 
for the exp-6 system though, actually, a careful examination of the data shows that the 
isostructural transitions, including the weaker ones, are many more. 3) In a wide pressure 
interval, the orthorhombic A20 lattice gives the most stable phase. There are actually five 
distinct A20 states at T = with abrupt transitions between them at specific pressures. 
These transitions are in the form of unequal compressions along the b and c axes, also 
accompanied by an adjustment of the internal parameter y. 4) The change in slope which 
is manifest in all curves at P ~ 52000 is due to a sudden change of the nearest-neighbour 
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distance of the optimal fee crystal from roughly 0.406 down to 0.37381 = Om (a similar jump 
occurs for bee at the same pressure). This effect is rather specific to the exp-6 system and 
ultimately related to the abrupt change of the potential profile at a M . 

Table 3 reports the T = phase diagram of the Yoshida-Kamakura potential. A com- 
parison with the exp-6 case reveals a number of similarities as for the trend of coordination 
number with pressure and for the order of appearance of the phases. Hence, some of the 
considerations made for the exp-6 model also apply to this potential, in spite of the fact 
that the steepness of the Yoshida-Kamakura potential is a monotonous function of r (see 
Fig. 1). In Fig. 4, we show the discrete RDF of the T = phases for uyk- We see that, with 
the exception of the highest-pressure phases, the range of distances from a central particle 
that corresponds to the force "plateau" is void of neighbours, like as if particles tended to 
sit at the shortest distances available at same force strength. Something similar occurs for 
the exp-6 model, see Table 2, where no particles are found in the region of distances where 
the potential is concave. 

We now try to obtain a rough phase diagram for the Yoshida-Kamakura model from 
just the knowledge of its T = sector. The simplest way to do this is to lay down a 



(Lennard- Jones-Devonshire) cell theory 45] for the solid phases and to use the Lindemann 



criterion 



46j for locating the melting transition. In the cell theory, a crystal partition 
function of effectively independent particles is written down where any given particle, which 
can be found anywhere in its own Wigner-Seitz cell (WSC), is acted upon by the force exerted 



by the other N — 1 particles, placed at their equilibrium lattice positions. In practice, 
where 



4 / d V- - / d%exp|-^0(r,)/(A; B r)l , (4.11) 

^WSCi ^WSCjv ^ j J 
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(we denote by capital letters the positions of lattice sites in the perfect crystal). Within this 
theory, the mean square displacement of a particle is given by 

/ Ar2 \ = / wsc d 3 rr 2 exp{-($(r)-$(0))/(/c i? r)} 

^ ' / wgc d3rexp{-(<&(r)-<&(0))/(fcBT)} ' K ' ] 
for $(r) = Ylj^i + r — ^ we denote r^^ the nearest-neighbour distance, the 
Lindemann fraction is defined as 

C = ^ '-, (4.14) 

12 



which is an increasing function of temperature. The Lindemann rule then states that melting 
occurs at the temperature T m where C reaches a critical value C c specific of the given lattice. 
Experience with other models shows that C c = 0.15-^-0.16 for fee and 0.18^-0.19 for bec 47] 
while, to our knowledge at least, there is no general consensus in the literature as to the 
value of C c for other lattices (we tentatively assume C c = 0.19). 

Fig. 5 shows the phase diagram of the Yoshida-Kamakura potential as mapped out in the 
way just explained. Aside from the simplicity of the theory used, in plotting this figure we 
make two further assumptions: 1) We discard the possibility that some other phases, not 
stable at T = 0, might be promoted entropically for T > 0; 2) we exclude that a given 
phase might invade, for T > 0, the density intervals where the adjacent phases are stable 
for T = 0. Looking at Fig. 5, we see at least one region of reentrant melting between bcc-I 
and sh-I, but probably there are others. This is a curious finding in view of the absence 
of two clearly-defined repulsive length scales for this potential (by contrast, reentrant-fluid 
behaviour is well-documented for the exp-6 model 9] and the Fomin potential 24| . both 
strongly-soft potentials). Since we do not calculate the phonon spectrum of the solids, we 
cannot exclude that some of them are actually mechanically unstable for T > 0. However, 
indirect clue to mechanical stability of a crystal phase can be the positive value of the elastic 
constant K associated with the one-particle potential $(r). This quantity can be extracted 
from the 0(r 2 ) term in $(r), via a spherical average 23j. We have checked that, at least 
within the range of stability of each phase in Table 3, the value of K is always positive. 
Moreover, if we estimate the mean square displacement of a particle from the harmonic 
approximation for $(r), the resulting phase diagram comes out not too different from that 
of Fig. 5. 

The T = phase diagram of the Fomin potential is outlined in Table 4. We see the 
same trend of z with pressure as observed in the previous cases. However, we find only one 
non-Bravais phase in this case (/3Sn) which turns out to be almost degenerate with the bet 
phase throughout the whole pressure range from 3.5 to 5.5. All in all, the phase behaviour 
of unbounded soft-core potentials has some recurrent features (e.g. rich solid polymorphism, 
low-coordinated non-Bravais crystal phases, and reentrant melting) which are also found in 
the phase diagram of many simple elements under extreme conditions. This suggests that 
the effective two-body (adiabatic) potential of these substances is a soft-core potential. 

Finally, we analyse the T = phase diagram of the Hertz potential, Eq. (12.41) . as sum- 
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marized in Table 5. A rather comp 
more complex than reported in 



ex behaviour shows up for this bounded potential, much 
251 ] . especially if compared with the simple phase portrait 
of the GCM. In fact, the phases listed in Table 5 are only those stable for P < 400 since 
apparently the sequence of T = phase transitions never comes to an end (as observed in 



Ref. 25j, there is no room for clustering in the Hertz model). This unique behaviour, never 
documented before, is the effect of a complex interplay between energy and volume consid- 
erations in the minimization of enthalpy as a function of pressure. Probably, this behaviour 
is related to the absence of a force plateau at r = which obliges the system to continuously 
setting right the positions of the neighbouring particles. 



V. CONCLUSIONS 



In recent years, increasing attention has been devoted to soft-matter systems as examples 

of anomalous thermodynamic behaviour, both in and out of equilibrium. This field of 

I I 

research is very much alive, with many points of contact with high-pressure physics [9|. We 
have been focused here on the anomalously rich solid polymorphism of systems of softly- 
repulsive particles, by studying the zero-temperature phase diagram of a number of model 
pair potentials with various forms of core softening. We found elements of complexity 
that are simply unknown to "normal" fluids interacting through a Lennard- Jones type of 
potential, with many low-coordinated non-Bravais lattices providing the structure of stable 
phases at intermediate pressures. In a near future, with the advent of new techniques to 
functionalize the colloidal surface, one can expect to obtain colloidal systems whose soft- 
core potential yields spontaneous assembling into similar exotic lattices. The relation of 
solid polymorphism to other kinds of thermodynamic oddities, such as reentrant melting 
and water-like anomalies, is the subject of work in progress. 

Acknowledgements 

We acknowledge useful discussions with Ezio Bruno and Sandro Scandolo. We also thank 
an anonymous Referee of our previous paper [9J for driving our attention to the papers 



14 



21 



22 



23]. 



[*] E-mail: Santi.Prestipino@unime.it 

[f] E-mail: saija@me.cnr.it 

[fl E-mail: malescio@unime.it 

[1] C. N. Likos, Phys. Rep. 348, 267 (2001). 

[2] G. Malescio, J. Phys.: Condens. Matter 19, 073101 (2007). 

[3] M. Watzlawek, C. N. Likos, and H. Lowen, Phys. Rev. Lett. 82, 5289 (1999). 

[4] C. N. Likos, N. Hoffmann, H. Lowen, and A. A. Louis, J. Phys.: Condens. Matter 14, 7681 
(2002). 

[5] A. A. Louis, P. G. Bolhuis, J. P. Hansen, and E. J. Meijer, Phys. Rev. Lett. 85, 2522 (2000). 

[6] F. H. Stillinger, J. Chem. Phys. 65, 3968 (1976). 

[7] A. Lang, C. N. Likos, M. Watzlawek, and H. Lowen, J. Phys.: Condens. Matter 12, 5087 
(2000). 

[8] S. Prestipino, F. Saija, and P. V. Giaquinta, Phys. Rev. E71, 050102(R) (2005). 

[9] G. Malescio, F. Saija, and S. Prestipino, J. Chem. Phys. 129, 241101 (2008). 

[10] A. B. Migdal, Sov. Phys. JETP 7, 996 (1958). 

[11] R. A. Buckingham, Proc. R. Soc. London, Ser. A 168, 264 (1938). 

[12] M. Ross and A. K. McMahan, Phys. Rev. B 21, 1658 (1980). 

[13] F. Saija and S. Prestipino, Phys. Rev. B 72, 024113 (2005). 

[14] M. Ross, R. Boehler, and P. Soderlind, Phys. Rev. Lett. 95, 257801 (2005). 

[15] See e.g. M. I. McMahon and R. J. Nelmes, Chem. Soc. Rev. 35, 943 (2006). 

[16] D. Gottwald, C. N. Likos, G. Kahl, and H. Lowen, Phys. Rev. Lett. 92, 068301 (2004). 

[17] D. Gottwald, G. Kahl, and C. N. Likos, J. Chem. Phys. 122, 204503 (2005). 

[18] J. Fornleitner and G. Kahl, Europhys. Lett. 82, 18001 (2008). 

[19] G. J. Pauschenwein and G. Kahl, Soft Matter 4, 1396 (2008). 

[20] G. J. Pauschenwein and G. Kahl, J. Chem. Phys. 129, 174107 (2008). 

[21] T. Yoshida and S. Kamakura, Prog. Theor. Phys. 47, 1801 (1972). 

[22] S. Kamakura and T. Yoshida, Prog. Theor. Phys. 48, 2110 (1972). 

[23] T. Yoshida and S. Kamakura, Prog. Theor. Phys. 56, 330 (1976). 

15 



[24] Yu. D. Fomin, N. V. Gribova, V. N. Ryzhov, S. M. Stishov, and D. Frenkel, J. Chem. Phys. 

129, 064512 (2008). 
[25] J. C. Pamies, A. Cacciuto, and D. Frenkel, cond-mat/0811 . 2227 
[26] S. Prestipino, F. Saija, and P. V. Giaquinta, J. Chem. Phys. 123, 144110 (2005). 
[27] See e.g. D. Frenkel and B. Smit, Understanding molecular simulation (Academic, New York, 

2001). 

[28] J. Behler, R. Martonak, D. Donadio, and M. Parrinello, Phys. Rev. Lett. 100, 185501 (2008). 
[29] T. Ishikawa, A. Ichikawa, H. Nagara, M. Geshi, K. Kusakabe, and N. Suzuki, Phys. Rev. B 

77, 020101(R) (2008). 
[30] A. R. Oganov and C. W. Glass, J. Chem. Phys. 124, 244704 (2006). 

[31] Y. Yao, J. S. Tse, Z. Song, D. D. Klug, J. Sun, and Y. Le Page, Phys. Rev. B 78, 054506 
(2008). 

[32] N. L. Abraham and M. I. J. Probert, Phys. Rev. B 77, 134117 (2008). 

[33] Useful information about a wealth of crystal lattices can be found at 

http : //cst-www . nrl . navy . mil/lattice/ index . html 
[34] M. I. McMahon and R. J. Nelmes, Phys. Rev. B 47, 8337 (1993). 

[35] U. Schwarz, K. Takemura, M. Hanfland, and K. Syassen, Phys. Rev. Lett. 81, 2711 (1998). 
[36] A. W. Overhauser, Phys. Rev. Lett. 53, 64 (1984). 

[37] M. Hanfland, K. Syassen, N. E. Christensen, and D. L. Novikov, Nature 408, 174 (2000). 
[38] M. I. McMahon, E. Gregoryanz, L. F. Lundegaard, I. Loa, C. Guillaume, R. J. Nelmes, A. K. 

Kleppe, M. Amboage, H. Wilhelm, and A. P. Jephcoat, Proc. Natl. Acad. Sci. U.S.A. 104, 

17297 (2007). 

[39] E. Gregoryanz, L. F. Lundegaard, M. I. McMahon, C. Guillaume, R. J. Nelmes, and M. 

Mezouar, Science 320, 1054 (2008). 
[40] J. B. Neaton and N. W. Ashcroft, Nature 400, 117 (1999). 

[41] N. E. Christensen and D. L. Novikov, J. Phys.: Condens. Matter 14, 10879 (2002). 
[42] H. Fujihisa, Y. Nakamoto, K. Shimizu, T. Yabuuchi, and Y. Gotoh, Phys. Rev. Lett. 101, 
095503 (2008). 

[43] A. N. Jackson, A. D. Bruce, and G. J. Ackland, Phys. Rev. E 65, 036710 (2002). 
[44] H. Eshet, F. Bruneval, and M. Parrinello, J. Chem. Phys. 129, 026101 (2008). 
[45] J. E. Lennard-Jones and A. F. Devonshire, Proc. Roy. Soc. A 163, 53 (1937). 

16 



[46] F. A. Lindemann, Z. Phys. 11, 609 (1910). 

[47] F. Saija, S. Prestipino, and P. V. Giaquinta, J. Chem. Phys. 124, 244504 (2006). 



17 




FIG. 1: Yoshida-Kamakura potential for a = 2.1 (solid line, in units of e). The figure reports 
also the force (dashed line, in units of e/cr) and the second derivative of the potential (dotted 
line, in units of e/cr 2 ). Note that the force and its steepness have been divided by two and four, 
respectively, in order to fit into the picture. 
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TABLE I: The non-Bravais lattices that have been considered in our search for stable phases at 
zero temperature. For each structure (column 1), the parameter (s) which are needed for a complete 
specification of the lattice are reported in column 2. The symbols in this column have the same 
meaning as in the web site of Ref. 33|], except for Imma 34j] and oC16-Cs [35]. The number of 
inequivalent sites (NIS, see main text) of the given lattice is also indicated in column 3. For Bravais 
lattices, the parameters are axial ratios (b/a,c/a) and angles (a, j3, 7), as usual. 



lattice 


parameters 


NIS (with weights, if different) 


hep 




1 


dhep 




2 


thep 




2 (1 x 1,2 x 2) 


9R 




2 (1 x 1,2 x 2) 


diamond 




1 


fluorite 




2 (1 x 1,2 x 2) 


A15 




2 (1 x 1,2 x 3) 


bccl2 




1 


/3Sn 


c/a 


1 


cI16-Li 


X 


1 


BC8 


X 


1 


graphite 


c/a 


2 


wurtzite 


(c/o)/V873,«/(3/8) 


1 


Imma 


c/a, A 


3 (1 x 1,2 x 1,3 x 2) 


BIO 


c/a, z 


1 


A7 


b/a, u 


1 


A8 


c/a, x 


1 


A13 


x 1 ,x 2 


2 (1 x 2,2 x 3) 


C6 


c/a, z 


2 (1 x 1,2 x 2) 


A20 


b/a,c/a,y 


1 


oC8-Ga 


(6/a)/1.69479, c/a, it/0.1549, u/0.0810 


1 


oC16-Cs 


(6/a)/0.594, (c/a)/0.590, x/0.2118, y/0.1781, z/0.328 


2 


ST12 


c/a, x x /0.0912, x 2 /0.1730, y 2 /0.3784, z 2 /0.2486 


2 (1 x 1,2 x 2) 
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TABLE II: Exp-6 potential for a = 11: Zero-temperature phase diagram. For each phase (column 
2), we show the pressure interval of stability (column 1), the related density interval (column 3), 
the first- and second- neighbour distances (column 4), the number of first- and second- neighbours 



(column 5), and the values of structure parameters (column 6). 



Pi-P 2 (10 3 eo-~ 3 ) 


phase 


P1-P2 {<j 3 ) 


n (0-) 

Ti (cr) 


Zl 
Z2 


parameters 


0-5.5 


fcc-I 


1.627-6.943 


0.58838-0.95435 

n oQonn 1 qac\pp 
0.o32U9-l. 34966 


12 

p 
b 


— 


5.6-11.8 


bcc 


6.999-9.267 


0.51947-0.57042 
59983-0 65867 


8 
p. 


— 


11.9-17.1 


hRl 


10.123-11.500 


0.37381 
0.53305-0.56628 


2 
6 


-0.227-=- -0.220 


17 on A 
1 1 .Z— ZU.4 


OUS-1 


ll.yoz— iz.o to 


0.37381 
0.52352-0.54667 


~ 3 
2 


0.85 -7- 0.91, 1.24 -7- 1.38, 
0.61^0.63, 2.06^2.24 


20.5-33.3 


/3Sn 


13.455-15.445 


0.37381 
0.52989-0.57802 


4 

4-5-12 


0.758 -7- 0.986 


33.4-36.9 


0C8-II 


16.464-16.722 


0.37381 
0.50274-0.50809 


~ 5 
2 


1.10 -7- 1.15, 1.81 -T- 1.85, 
1.02 -=- 1.04, 1.30 -7- 1.34 


37.0-45.3 


0C8-III 


18.561-18.651 


u.o 1 001 
D ^4^46-0 55677 


~ 7 



95 — 98 1 08 — 1 09 
1.10 -7- 1.12, 1.17 -7- 1.23 


45.4-47.6 


0C8-IV 


19.097 


0.37381, 0.51223 


~ 7, 2 


0.88, 1.11, 1.17, 1.33 


47.7-48.2 


A20-I 


19.642 


0.37381, 0.56016 


~ 8, 4 


1.726, 0.666, 0.167 


48.3-52.1 


A20-II 


19.905 


0.37381, 0.54115 


~ 8, 4 


1.522, 0.629, 0.179 


52.2-52.8 


A20-III 


20.074 


0.37381, 0.53498 


~ 8, 4 


1.468, 0.621, 0.183 


52.9-57.3 


A20-IV 


20.383 


0.37381, 0.52665 


~ 8, 4 


1.398, 0.610, 0.189 


57.4-60.9 


A20-V 


20.813 


0.37381, 0.51700 


~ 8, 4 


1.327, 0.601, 0.196 


61.0-68.3 


sh 


22.107 


0.37381, 0.52864 


8, 12 


1 


68.4-76.9 


wurtzite 


24.337 


0.37381, 0.52864 


10, 9 


2.225, 0.734 


77.0- 


fcc-II 


27.075 


0.37381, 0.52864 


12, 6 
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TABLE III: Yoshida-Kamakura potential for a = 2.1: Zero-temperature phase diagram. Notations 
are the same as in Table 2. 



P1-P2 (ea- 3 ) 


phase 


P1-P2 (^" 3 ) 




z 


parameters 


0-0.63 


fcc-I 


0-0.325 


1.63259- 


12 


— 


0.64-1.26 


bcc-I 


0.331-0.415 


1.46282-1.57736 


8 




1 97 9 90 


cVi T 

sn-i 


U.400— u.oy4 


u.yu/oo— u.yyoio 


_ 


U.OIO — U.040 


9 qn 9 
z.ou— Z.OD 


A 7 


r\ p.-\ A n RQ7 
U.014— U.Oo 1 


u.yo4io— u.yoy < 1 


q 



3.962 — 3.975, 0.181 — 0.184 


z.od— 4.» 1 


diamond 


n R7n n S/i k 

U.O 1 U— U.o4D 


u.yiouo— u.yoy < u 


A 

4 




A Q9 ^ /Ifi 
4.y/— 0.40 


qVi tt 
sn-11 


n ons n Q/i/i 
u.yuo— u.y44 


n Q798Q 

u.yoo ( 4— u.y / zot) 





l.ool — l.ooo 


5.47-12.07 


A20 


0.985-1.322 


0.90053-0.97273 


~ 8 


1.728 -7- 1.731, 0.626 0.646, 
0.167 


12.08-15.68 


hcp-I 


1.424-1.583 


0.96311-0.99770 


12 




15.69-52.75 


bcc-II 


1.593-2.618 


0.79168-0.93426 


8 




52.76-138.28 


fcc-II 


2.629-4.024 


0.70570-0.81328 


12 




138.29-365.65 


hcp-II 


4.026-6.075 


0.61516-0.70558 


12 




365.66- 


fcc-III 


6.077- 


-0.61509 


12 





TABLE IV: Fomin potential for ro = 1.35 a: Zero-temperature phase diagram. Phases with 5 
parameters were not included in our scrutiny. Notations are the same as in Table 2. 



P1-P2 (ecr- 3 ) 


phase 


P1-P2 (o-~ 3 ) 




z 


parameters 


0-3.12 


fcc-I 


0-0.468 


1.44574- 


12 




3.13-3.55 


bco 


0.549-0.559 


1.07840-1.08786 


2 


1.326 ~r~ 1.330, 2.134 -=-2.145 


3.56-4.62 


bet 


0.655-0.680 


1.08276-1.09938 


4 


2.298^-2.317 


4.63-5.54 


/3Sn 


0.690-0.708 


1.05191-1.06223 


4 


0.729 -j- 0.733 


5.55-8.91 


sc 


0.801-0.859 


1.05197-1.07677 


6 




8.92-14.88 


sh 


0.940-1.017 


1.03906-1.06239 


~ 8 


0.988 -7- 0.994 


14.89- 


fcc-II 


1.166- 


-1.06645 


12 
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TABLE V: Hertz potential: Zero-temperature phases up to P = 400. For P > 40, only phases 
with and 1 parameters were examined. Notations are the same as in Table 2. 



Pt-P 2 (ea- 3 ) 


phase 


PI-P2 (cr^ 3 ) 


n (a) 


z 


parameters 


0-0.49 


fcc-I 


0-2.211 


0.86161- 


12 




0.50-1.70 


bcc-I 


2.287-3.447 


0.72232-0.82817 


8 




1.71-2.45 


sh-I 


3.694-4.231 


0.57655-0.60419 


2 


0.838 -7- 0.840 


2.46-3.96 


ti n t 

C116-1 


4.353-5.249 


0.51263-0.54563 


3 


0.125 


3.97-5.73 


pbn-1 


5.485-6.386 


0.50503-0.52598 


4 


2.416 — 2.508 


5.74-7.20 


X. 4- T 

bct-1 


r* A r -7 11/'' 

6.456-7.116 


n -1 rn'7 A n /irT^i 

0.45974-0.46571 


2 


0.571 — 0.588 


7. 21-10. bl 


A OA T 

A20-1 


7.366-8.651 


0.48857-0.51514 


~ 6 


1.730 — 1.731, 0.719 — 0.728, 0.167 


10.62-11.10 


fluorite 


8.745-8.910 


0.47819-0.48118 


8x 1 
4x2 


— 


11.11-11.57 


fcc-II 


9.050-9.214 


0.53541-0.53863 


12 




11.58-14.63 


A20-II 


9.331-10.309 


0.49960-0.51378 


~ 6 


1.731, 1.448 -j- 1.475, 0.167 


14.64-28.14 


bcc-II 


10.450-13.895 


0.45386-0.49908 


8 




28.15-28.64 


psn-n 


14.093-14.202 


0.40956-0.41061 


4 


0.578 


28.65-38.10 


sh-II 


14.254-16.214 


0.38467-0.40155 


6 


0.894-^0.895 


38.11-54.07 


sc 


16.417-19.265 


0.37303-0.39346 


6 




54.08-58.89 


hRl-II 


19.431-20.218 


0.37193-0.37660 


6 


-0.074^-0.072 


58.90-75.02 


bct-II 


20.312-22.789 


0.32347-0.33647 


2 


0.617^-0.622 


75.03-102.78 


hRl-III 


22.951-26.643 


0.34178-0.36507 


6 


0.543 -7- 0.553 


102.79-141.62 


hep 


26.905-31.323 


0.35608-0.37459 


12 




141.63-148.08 


bct-III 


31.413-32.096 


0.34382-0.34639 


8 


0.916^0.923 


148.09-169.07 


bct-IV 


32.121-34.248 


0.33691-0.34488 


8 


1.116 1.155 


169.08-206.32 


bct-V 


34.328-37.754 


0.31071-0.32268 


4 


1.734 1.766 


206.33-242.37 


bccl2 


37.822-40.873 


0.31085-0.31900 


8 




242.38-307.75 


hRl-IV 


41.044-46.081 


0.28614-0.30195 


6 


0.130^0.171 


307.76-345.22 


hRl-V 


46.113-48.746 


0.27802-0.28422 


6 


0.098^0.110 


345.23-370.23 


/3Sn-III 


48.859-50.541 


0.27074-0.27429 


4 


3.967^-3.988 


370.24-400 


cI16-II 


50.669-52.602 


0.25402-0.25721 


3 


0.078 
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FIG. 2: (Color online) Effective inverse-power exponent n(r) for the potentials (|2.ip - (|2.4p . see 
Eq. (|2.5h : Exp-6 potential (solid line), Yoshida-Kamakura potential (dotted line), Fomin potential 
(dashed line), and Hertz potential (long-dashed line). For the exp-6 case, n(r) was plotted only 
for those r where ue{t) > 0. 
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FIG. 3: (Color online) Exp-6 potential for a = 11: T = chemical potential fi, plotted as a 
function of pressure P, for a number of crystal structures, choosing the fee lattice as reference 
(both and P are in reduced units; the chemical potentials of structures that are never stable 
are not shown). Besided fee, the stable phases are bec (squares and solid line), hRl (triangles and 
solid line), /3Sn (dots and dotted line), oC8 (solid line), A20 (crosses and solid line), sh (dashed 
line), and wurtzite (starred crosses and red dashed line). 
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fcc-I 




bcc-I 1 


sh-I . 1 1 . 


A7 . . 1 


diamond 

i 




sh-II , .11 


A20 

,i ii 


1, , 


hcp-I . 


1 


bcc-II | 




1 


1 
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1 


1 1 1 1 1 1 1 1 i I""" — r- T " ~T 1 



0.5 1 1.5 2 

r/ o 



FIG. 4: (Color online) Yoshida-Kamakura potential: Comparison of the discrete radial distribution 
function for the various zero-temperature phases (the highest-pressure hcp-II and fcc-III phases 
are not shown). At every distance, this function simply counts the number of neighbours at 
that distance. The height of each panel is 15. Bottom panel: the potential (solid line) and the 
corresponding force divided by two (dashed line). From top to bottom: fcc-I (P = 0.32), bcc-I 
(P = 0.95), sh-I (P = 1.79), A7 (P = 2.43), diamond (P = 3.73), sh-II (P = 5.19), A20 (P = 8.77), 
hcp-I (P = 13.88), bcc-II (P = 34), and fcc-II (P = 56). 
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FIG. 5: (Color online) Yoshida-Kamakura potential: Schematic phase diagram in the p-T plane up 
to per 3 = 3, as obtained by combining the cell theory for the T = solids with the Lindemann rule. 
Density gaps between the phases are two-phase coexistence regions. In the low-density region, the 
prediction of Eq. (|2.10p for the fee melting line (open dots) is compared with that based on the 
Lindemann criterion. 
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